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We analyze magnetic flux tubes at zero temperature in a superconductor that is coupled to a 
superfluid via both density and gradient ( "entrainment" ) interactions. The example we have in 
mind is high-density nuclear matter, which is a proton superconductor and a neutron superfluid, 
but our treatment is general and simple, modeling the interactions as a Ginzburg-Landau effective 
theory with four-fermion couplings, including only s-wave pairing. We numerically solve the field 
equations for flux tubes with an arbitrary number of flux quanta, and compare their energies. 
This allows us to map the type-I/type-II transition in the superconductor, which occurs at the 
conventional « = A/£ = l/s/2 if the condensates are uncoupled. We find that a density coupling 
between the condensates raises the critical k and, for a sufficiently high neutron density, resolves 
the type-I/type-II transition line into an infinite number of bands corresponding to "type-II(n)" 
phases, in which n, the number of quanta in the favored flux tube, steps from 1 to infinity. For 
lower neutron density, the coupling creates spinodal regions around the type-I/type-II boundary, in 
CO \ which metastable flux configurations are possible. We find that a gradient coupling between the 

condensates lowers the critical k and creates spinodal regions. These exotic phenomena may not 
occur in nuclear matter, which is thought to be deep in the type-II region, but might be observed 
in condensed matter systems. 

b 

' I. INTRODUCTION 

a/ 

Superconductivity and superfluidity are well-studied phenomena, known to occur in many physical systems, from 
s ! ' cold metals and cold atomic gases to nuclear matter and quark matter. In this paper we investigate a system that 
has both a charged condensate, leading to superconductivity, and a neutral condensate, leading to superfluidity. We 
focus on the magnetic flux tubes that are associated with the superconducting condensate, and study how they are 
modified by the presence of the superfluid, assuming that the two condensates can interact with each other via density 
and gradient ("entrainment") interactions. 

An example of this type of system is nuclear matter, which at sufficiently high density undergoes Cooper pairing of 
both neutrons and protons. We will present our calculations in this context, referring to the charged condensate as 
the "proton condensate" and the neutral one as the "neutron condensate" , and choosing values appropriate to nuclear 
matter for our parameters when presenting numerical results. In fact, the questions that we study in this paper were 
originally raised in investigations of the nature of the proton superconductivity in the nuclear matter in a neutron 
star. Although it is generally believed that the protons form a type-II superconductor [1], there is evidence from long 
neutron star precession periods that seems to favor type-I superconductivity [2] (for contrary views see [3, 4]). This 
led Buckley et. al. [5] to suggest that, if the density interaction between the magnitudes of the neutron and proton 
Cooper pair condensates is extremely strong, nuclear matter would be a type I superconductor even if its penetration 
depth A and coherence length £ obey the conventional condition A/£ > l/v2 for type II superconductivity. We have 
argued elsewhere that the assumption of a strong coupling between the proton and neutron condensates is wrong for 
neutron star matter [6]. However, Buckley et. al. were correct in making the point that a superconductor will be 
affected by interactions with a co-existing superfluid. 

In this paper we study the type I versus type II nature of a (proton) superconductor coupled to a (neutron) 
superfluid, using an effective theory for the protons and neutrons that contains four-fermion interaction terms which 
lead to s-wave pairing. We do not include highcr-angular-momcntum pairing, although that would be needed for a 
more realistic analysis of high-density nuclear matter. Our analysis extends that of Ref. [5] in the following ways: 
(a) Our model, like that of Ref. [5], contains a coupling a np between the magnitudes of the neutron and proton 
condensates, and self-couplings a nn and a pp , but we survey the whole range of values of a np , from zero to of order 
a pp ; (b) we also include "entrainment" interactions between the gradients of the proton and neutron condensates; 
(c) we use a simpler and more direct method to study the type-I/type-II phase boundary, using the energetics of 
flux tube coalescence/fission: we calculate the energy of flux tubes with a wide range of magnetic fluxes, from one 
quantum to several hundred quanta, and find which one has the lowest energy per unit flux. As we will see, this 
has the additional benefit of allowing us to find exotic stable multi-quantum flux tubes, such as have been found in 
systems of two coupled superconductors [7]. However, as we discuss below, our analysis is not sensitive to minima in 
the interaction energy at finite separation between flux tubes. 

Our analysis is entirely at zero temperature. This is a good approximation for neutron star matter near nuclear 
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saturation density, where the critical temperatures for the superfiuid and superconductor are of order MeV [8-10]. 
The temperature of a compact star drops below this value within minutes of its formation in a supernova, and is at or 
below the keV range after the first 1000 years [11]. When we discuss type-I versus type-II behavior we are referring 
to the response of the system to a magnetic field at the lower critical value, at T = 0. 

As far as we know, there has been no previous work on how a flux tube in a superconductor is affected by a 
gradient coupling to a co-existing superfiuid. However, there has been work on possible knot solitons [12], vortices 
in the SO(5) model of high-temperature superconductivity [13], and on the complementary situation, a superfiuid 
vortex with gradient coupling to a co-existing superconductor. There the coupling leads to the "entrainmcnt" or 
Andrccv-Bashkin effect [14] whereby the proton condensate is dragged along with the neutron condensate, producing 
a non-zero proton current around the vortex, dressing it with some magnetic flux [15]. It is interesting to note that 
this flux is not a multiple of the flux quantum for proton flux tubes. This is possible because of the difference between 
the energetics of a neutron vortex and a proton flux tube. The flux tube has energy density localized to the vicinity 
of its core. Far from the core the energy density must vanish, which means the proton field must change in phase 
by a multiple of 27r, and the the vector potential must cancel the resultant gradient, leading to a quantized magnetic 
flux. A neutron vortex, by contrast, has gradient energy that is not localized to the vicinity of the vortex, and the 
total energy per unit length diverges in the infinite volume limit. The vector potential is therefore not constrained to 
cancel any gradient in the proton field, and takes on a value that minimizes the overall energy, with no quantization 
condition on the resulting magnetic flux. 

Returning to the situation that we study, a proton flux tube in a neutron superfiuid background, we do not expect 
a similar behavior. This is because the proton flux tube's energy density is localized around its core, giving it (unlike 
the neutron vortex) a finite energy per unit length. If the neutron condensate were entrained, and developed non-zero 
circulation around the flux tube, it would acquire a non-localized energy density, leading to an infinite energy per unit 
length for the flux tube, which is clearly energetically disfavored. We will see below that the effect of gradient couplings 
on the proton superconductor is more subtle: it leads to mctastablc regions near the type-I/type-II boundary. 



II. STABILITY OF FLUX TUBES 



Our aim is to explore the response of the proton superconductor to an applied critical magnetic field at zero 
temperature. We will therefore construct a phase diagram in the space of the coupling constants of the Ginzburg- 
Landau effective theory. We would like to be able to specify when it is of type II (at the lower critical magnetic 
field, flux tubes appear, and remain separate, i.e they repel) and when it is of type I (at the critical magnetic field, 
macroscopic normal regions appear, i.e. the flux tubes attract and coalesce). The simplest way to do this is to 
calculate the energy per unit length E n of a flux tube containing n flux quanta. The same approach has been used 
for vortices in the SO(5) model [16]. It is convenient to work in terms of the energy per flux quantum, 

B n = ^-E 1 . (1) 
n 

When B n is negative the n-quantum flux tube is stable against fission into many single quantum flux tubes, and it is 
energetically favorable for n single quantum flux tubes to coalesce into one n-quantum flux tube. When B n is positive 
the n-quantum flux tube is unstable against fission, and coalescence is energetically disfavored. If one calculates B n 
for all n then the energetically favored value of n is the one that minimizes B n . 

In a traditional type I superconductor, small flux tubes attract each other and amalgamate into large ones and 
ultimately into macroscopic normal regions, so we would expect to find B n < with its value dropping monotonically 
as n rises. In a type II superconductor we would expect B n > 0, with its value rising monotonically with n. Our 
calculations confirm these results for a single superconductor, but we will see that B n shows more complicated behavior 
when the superconductor feels interaction with a co-existing superfiuid. 

Calculations of B n are straightforward because they always occur in a cylindrically symmetric geometry, so the 
problem is one-dimensional. For a more detailed understanding of flux tube interactions, one would have to consider 
two single-quantum flux tubes a distance d apart. Their total energy is U(d), where U(0) — E 2 and ?7(oo) = 2Ei, so 
B 2 = i(L/(0)-[/(oo)). As expected, B 2 < means that the flux tubes have lower energy when they amalgamate, and 
B 2 > means that the flux tubes have lower energy when they separate. If U(d) is monotonic, we can conclude that 
flux tubes either coalesce (B 2 < 0) or repel to infinite separation (B 2 > 0), corresponding to type-I or type-II behavior 
respectively. However, if there is a minimum in U(d) at some favored intermediate separation d = d* then irrespective 
of the sign of B n , one has a new variety of type II superconductor with some favored Abrikosov lattice spacing d*. 
Such behavior has been found to arise from a (f> 6 term [17] and in the case of two charged condensates [7]. Calculating 
U{d) in the current context is an interesting but demanding problem which we leave for future work. In this paper 
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we assume that U(d) is monotonic, so to analyze the attractiveness/repulsiveness of the flux tube interactions it is 
sufficient to calculate B n , or equivalcntly E n /n. 



III. FLUX TUBES IN THE GINZBURG-LANDAU MODEL 



A. Ginzburg-Landau model 

We start by writing down the zero-temperature Ginzburg-Landau effective theory of proton and neutron condensates 
in the presence of a magnetic field [5, 18]. We denote the proton condensate field by 4> p , the neutron condensate field 
by <j) n , and the magnetic vector potential by A. The free energy density is 

^ = ^-(l(V - ^A)<M 2 + IV0J 2 ) + |V ^ + U ent (<f> p) <f> n ) + V{\j> p \\ \4>n\ 2 ) (2) 
2m c he 8tt 

where m c is twice the nucleon mass, q is twice the proton charge, U en t is the cntrainment free energy density (see 
[18]) 



U e 



2m c 2(4>p)((j) n ) 



(V - ^A)0 P • V0„j + <^„ ^(V - ^A)0 P • V# 

(v + ^A)0; • vc) + Mn (( v + y c a K ■ v<£„ ) 



and 



V(\^ P \ 2 , \M 2 ) = -VM 2 - Vn\K\ 2 + ^f\M 4 + ^Kl 4 + aM 2 ^ 2 (4) 

a is a parameter characterizing the strength of the gradient coupling, \x v and fx n are the chemical potentials of the 
proton and neutron condensate excitations, and a pp , a nn , and a pn are the GL quartic couplings. 

In zero magnetic field, the condensates would have position-independent bulk densities (4> v ) 2 and (<f> n ) 2 obtained 
by minimizing the free energy. This allows us to eliminate the chemical potentials fj, p , p n by writing 

p p = a pp (<j) p } 2 + a pn (<j) n } 2 

Mn = a nn ((/) n } 2 + a pn {(j) p ) 2 (5) 

so up to constants involving (<fi p ) and (4> n ), the potential V can be expressed in terms of the deviations of the 
condensate fields from their bulk values: 

t/(I0 p i 2 , |0„i 2 ) = y (W 2 r & P ) 2 )l+*r(\M 2 - <<M 2 ) 2 (6) 



-apn {\4> P \ 2 - {<t>p) 2 ) (!</>« 



12 _ I A \2^ 



In a neutron star, electrical neutrality keeps the proton fraction small, in the 5% to 10% range [19, 20]; we will 
take (4> p ) 2 1 '((fin) 2 ~ 0.05. As we now argue, a typical value for the entrainment coupling is a ~ 10" 1 . We first relate 
our formalism to the hydrodynamic limit of the free energy, following [18]. We focus on the phases of the fields, 
4>p = {4>p) gx p(*Xp) an d 4>n = (4>n) cxp(ixn), and assume the fields have constant magnitude, and their phases have 
gradients 

ft ^ 2e , h „ 

Vx P A , v„ = Vx„ ■ (7) 



2m p rripC 2m 

The free energy density (2) then reduces to the hydrodynamic form 

F = \ P pp v 2 p + l -p n ^ 2 n + p""v p ■ v„ + V + ^ , (8) 

where the symmetric matrix p of supcrfluid densities has elements 

p pp = 2m p ((f> p ) 2 » m c (0 p ) 2 , p"" = 2m„(0„) 2 » m c (0„) 2 , P pn = -2m n a(4> p )(cb n ) . (9) 

Our entrainment parameter a is therefore related to the parameter e of Ref. [21-23] by a = e((j> n ) / (<f> p ). Since e is of 
order 0.03, and (</>„) 2 / '(4> p ) 2 ~ 20, we expect a ~ 10 _1 . This is consistent with the estimate p p ™ ss —^p pp used by [18]. 
In terms of the Andreev-Bashkin parametrization [14], p\% = ~p pn , p\ = p pp + p pn , pi = p nn + p pn , so p\j pn ~ 1, 
P2I p\2 ~ 40. All the interactions in (8), including the entrainment, have their ultimate origin in the strong interaction 
between the nucleons, which is isospin symmetric, and hence does not distinguish protons from neutrons. 
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B. Flux tube solutions 

To study a flux tube containing n flux quanta, we assume a cylindrically symmetric field configuration in which the 
proton condensate field winds (in a covariantly constant way) around the z-axis with a net phase 27m, 



<t> P = (<P P ) f(r)e l 
K = (<f>n) g(r) 
_ nhca(r)~ 



(10) 
(11) 

(12) 



q r 

We have defined <j> n as a real field, because, as noted above, any net phase change in the neutron condensate when it 
circles the flux tube would cost an infinite energy per unit length. Inserting the ansatz in (2) we obtain 



T = 



& r,, , n 2 / 2 (l-a) 2 



2m r 



+ (cj> n ) 2 (g') 2 -2*(</> p )(<f> n )f-g-f'-g' 



n 2 h 2 C 2 {a') 2 CLppjfip) 4 ( f 2_-,\2, OnJA^r r 2 _ -, y 
8Tiq 2 r 2 2 U ' + 2 [9 > 

+ a pn ^ p ) 2 (<f> n ) 2 (f 2 ~l) ( 5 2 -l) 



(13) 



Generating the Euler-Lagrange equations using the standard procedure, we obtain a set of coupled differential 
equations for /, g and a: 



2m c a pp (4> p ) 2 



f „ + f_ » 2 (l~a) 2 / 



[f-g(g" + g -)+f{g') 



2vn c (x pp ( 



m c c 



Airq 2 



r 



n a 
a 

r 



/(/ 2 -l) + 



-f(9 2 1) 



[ /" + 7 ) 



^ 2 -i) + ^s(/ 2 -i) 



'!>!> 



(l-«)/ 2 



(14) 



At this point we recall the definition of the Ginzburg-Landau parameter k = A/£, where the London penetration 
depth A and superconducting coherence length £ are (see [24]) 



A = 



m c c A 



m c cr 



Anq 2 (<j) p ) 2 



\QnhcaEM{<j)p) 2 



n 2 



2m c a pp {4> p ) 2 



(15) 



To further simplify the equations, we then change variables to a dimensionless radial coordinate f — r/£, obtaining 



f n\\-a) 2 ! i 
J + — ^ o- 

r r i 



[g" + j)+.f{g'f 
\f"+ f j)+g{f? 



= /(/ 2 -i) + 



A'pn Yrn/ 



:f(g 2 i) 



^ 2 -i) + ^ 3 (/ 2 -i) 



a"-% = -4,(1 -a)/ 2 
The free energy per unit length of the flux tube, in terms of the variable f, is 



2ira z 



i 2 2 
+71 K 



(a') 2 , 1 
f 2 2 



+ ^(/ 2 -l) + 



2 la 



2 a 



PP vrpl 



~(9 2 -l) 2 ' 



*pp Wpl 



(f 1) (9 2 1) 



(16) 



(17) 
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FIG. 1: (Color online) Profile of flux tube with n = 1 units of flux (left) and n = 100 units of flux (right) showing the effect 
of density coupling j3 between neutron and proton condensates. The plot shows the deviation Sp of the condensates from their 
vacuum values (18). With no coupling between the condensates (/3 = a = 0), the neutrons are undisturbed (Sp n = 0). With a 
non-zero density coupling f3, the neutron condensate (broken lines) is significantly perturbed by the flux tube. Note that the 
neutron Sp n 's are multiplied by 10 (not by 100 as in Fig. 2) to make them visible. The other parameters are k = 3.0, a = 0.0, 
and <0 p ) 2 /<0n) 2 - 05. 



In addition to the system of equations, we require boundary conditions on the fields at the origin and at oo. Far 
from the flux tube core, the fields will go to their uniform condensate value, so /(oo) = g(oo) = a(oo) = 1. Near the 
origin, f(r) oc r n , a(r) oc r 2 and g(r) is a constant. Therefore we have the conditions /(0) = 0, a(0) = and g'(0) = 0. 
To obtain the energy of a flux tube we numerically solve the ODE system for the neutron and proton condensate and 
magnetic potential profile functions, then calculate the free energy of the system by inserting the results into (17) and 
integrating. 

The system has five independent parameters: a pp , a nn /a pp , a pn /a pp , a, and (4> n ) / {<fip} ■ In neutral nuclear matter, 
the density of protons (neutrons) is proportional to (4> p ) 2 ((</>n) 2 ), and the proton density is approximately 5% of the 
total baryon number density [18], so we set (4>p) 2 / (4>n) 2 = -05 in most of our analysis. Following [5, 6] we set a nn = a pp , 
and use (15) to exchange the parameter a pp for k, which is the conventional parameter used in condensed matter 
studies of superconductivity. Our reduced set of parameters is therefore k, the proton-neutron gradient coupling er, 
and the proton-neutron density coupling (3 = a pn /a pp . We also study some effects of varying (4> P ) 2 / (0n) 2 - 



IV. NUMERICAL RESULTS 



A. Flux tube solutions 



For given values of (4> p ) 2 / (4>n) 2 , k, the proton-neutron gradient coupling a, and the proton-neutron amplitude 
coupling (3 = a pn /a pp we numerically solved the equations of motion (16) giving the field profiles for flux tubes with 
various numbers n of flux quanta. We obtained the solutions using a finite-element relaxation method, which is much 
less sensitive to initial conditions than the traditional "shooting" method, and better suited to repeatedly solving the 
equations for different sets of parameters. Next, we insert the solution for each profile into our expression for the free 
energy (17) and numerically integrate it to obtain a value for E n . 

To estimate the numerical errors in our results, we varied the convergence criterion in the finite-element relaxation 
calculation, the spacing of the radial grid of points, and the radius out to which the grid extended. We found that 
the resultant variation in E n /n was of order 10 -6 , so numerical errors are invisible on the scale of the plots shown in 
Fig. 3. 

Having obtained E n we then plot the series B n to determine whether the system is type I or type II for the chosen 
point in parameter space. In this way we find the points in parameter space where the system changes from a type 
I state to a type II state. Taking various slices through the parameter space, we can generate phase diagrams that 



G 



Profiles with nonzero gradient coupling 



i 



Single flux quantum 



5p p 

-- 5p n x 100(0 = 0.5) 
■-■ 5p n x 100(o = -0.5) 
i , i 



10 20 
radial distance r = r/^ 



30 



100 flux quanta 




5p p 

- - 5p n x 100 (o = 0.5) _ 
■ - ■ 5p„ x 100 (o = -0.5) - 



10 20 
radial distance f = r/i; 



30 



FIG. 2: (Color online) Profile of flux tube with n = 1 units of flux (left) and n = 100 units of flux (right) showing the effect 
of gradient coupling a between neutrons and protons. The plot shows the deviation 8p of the condensates from their vacuum 
values (18). With no coupling between the condensates (/3 = a — 0), the neutrons are undisturbed (Sp n = 0). With a non-zero 
gradient coupling a, the neutron condensate (broken lines) is slightly perturbed by the flux tube. Note that the neutron <5p n 's 
are multiplied by 100 to make them visible. The other parameters are k = 3.0, (3 = 0.0, and (0 P ) 2 /<0„) 2 =.O5. 



show the boundary curves between the various phases. 

Figs. 1 and 2 each show a profile for a flux tube with a single flux quantum n — 1 on the left, and a profile for a 
flux tube with 100 flux quanta on the right. Fig. 1 shows the effect of non-zero density coupling j3 and Fig. 2 shows 
the effect of non-zero gradient coupling a. We have plotted the normalized difference in density of the pair fields from 
their condensate values, 

_ ^(f)-(0 p ) 2 
S Pp\ r ) = 72 \ 2 = f ( r ) ~ 1 



6Pn(r) - €( t- { f n)2 =9 2 (r)-l (18) 



1. No coupling to neutrons 



We do not show a plot of the flux tube profile for a simple superconductor, since this is well known: in a core region 
whose area rises as the number of flux quanta n, the proton condensate is suppressed; in a wall region the condensate 
returns to its vacuum value. At the Bogomolnyi point [25], k = l/y/2, the energy per flux quantum is independent of 
n [26], but on cither side of this value there are area and perimeter contributions to the energy [27], so for k close to 
l/y/2 we expect the energy of a flux tube in a simple superconductor to have the following dependence on n, 

E { n sc \n) = nE Bog +8kM(u- ci \AI + ci + •• ■) . (19) 

This is an expansion around n = oo, but our numerical results will show that it works down to n = 1. Wc define 
0k = k — l/y/2. E Bog is the energy per unit flux at 5k = 0. By convention we take the parameter M, which has 
dimensions of energy, to be positive. The value of ci is then positive, ensuring that for Sk > 0, n = oo is disfavored 
(type-II), and for 6k < 0, n = oo is favored (type-I). We will see this behavior in our numerical results (Sec. IV B 1 
and upper left plot of Fig. 3). 
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2. Density coupling to neutrons 

For positive /3, which corresponds to positive a np , equations (2) and (4) indicate that there is a repulsion between 
the neutron and proton condensates, so in the center of the flux tube, where the proton condensate is suppressed, the 
neutron condensate will be enhanced. That is exactly what we see in Fig. 1, where the dashed curve, showing the 
perturbation to the neutron density p n , rises inside the flux tube. For negative (3 there is attraction between the two 
condensates, and the neutron condensate is suppressed inside the flux tube (dash-dotted line). We therefore expect 
that the leading correction due to the interaction will be proportional to the core area, i.c proportional to n. The 
energy of an n-quantum flux tube is then 

E n (n,0) ^E^(k) +M p (-n + bi y/n + &! + ■••) , (20) 

where En^ («) is the energy for an n-quantum flux tube in a pure superconductor, with no coupling to a superfluid 
(19). The leading correction is —Mpn, which should be negative and quadratic in (3 for small (3 (see Sec. IVA4), so 
the interaction energy parameter Mp is positive and proportional to f3 2 . The sub-leading term proportional to y/n 
arises from the energy cost of the gradient in p n at the edge of the flux tube, where it must return to its vacuum 
value, so we expect this term to be positive: bi > 0. We do not have an a priori expectation for the sign of the 
sub-sub-leading term b\. 

3. Gradient coupling to neutrons 

For positive a, we expect from (2) and (3) that the positive gradient in p p at the wall of the flux tube will induce a 
positive gradient in p n in the same range of radii, which lowers the energy of the system. This is exactly what we see 
in Fig. 2, where the dashed curve showing the perturbation to p n has a positive slope in the range of radii where the 
solid curve (p p ) has the largest positive slope. On either side of that region it has a negative slope, as it returns to 
its unperturbed value. For negative a the effect is reversed: the dash-dotted curve shows p n having a negative slope 
where p p has the largest positive slope. 

We therefore expect that in the presence of a gradient coupling, the correction to the energy of a flux tube has a 
dominant core-perimeter term proportional to \/n, 

E„(k, a) « E^ c \n) + M a (-si s/h~ + «! + •••)■ (21) 

The energy correction is negative and quadratic in a for small a (see Sec. IV A 4), so the interaction energy parameter 
M a is proportional to a 2 ; choosing it to be positive by convention requires si to be positive. We do not have an a 
priori expectation for the sign of si. 

4- Symmetry under change of sign of couplings 

It is clear from Figs. 1 and 2 that for couplings /3 and a of order 0.5 the modification of the field configuration due to 
the interaction between the condensates is extremely small, so it is reasonable to treat its effects perturbatively. (At 
the end of Sec. IV C we will discuss the limit of small neutron condensate, where the perturbativc approach becomes 
questionable.) 

When we evaluate the perturbative correction to the energy of the flux tube, there is no linear term in (3 and a. Such 
a term would arise from evaluating the (3 and a terms from the Hamiltonian in the unperturbed field configuration. 
But in that configuration the neutron condensate sits at its vacuum value, so both terms evaluate to zero (g = 1, 
g' = in (17)). 

We therefore expect the change in the energy of the flux tube to be quadratic in the couplings (3 and a. Firstly, 
this correction must be negative. This is a well-known result from perturbation theory: the second-order correction 
arises from the change in the configuration in response to the perturbation, which only occurs because it is driven by 
a resultant lowering of the energy. Secondly, the change in the energy will in general contain (3 2 , a 2 , and (3a terms. 
This means it will be even in (3 when a = and even in a when (3 = 0, so we expect Mp oc {3 2 and M a oc a 2 in 
Eqs. (20) and (21). 

However, if both (3 and a are nonzero, then the (3a terms spoil the symmetry of the energy under negation of the 
couplings. This is clear from Figs. 1 and 2. For example, suppose that as well as non-zero j3 we have a very small 
non-zero a. Now consider sending (3 — > — /3. From Fig. 1 we see that this changes the sign of the slope of p n in the 
wall region where p p has positive slope. If a is nonzero then these two configurations will have different energies, since 
the gradient of p n is then coupled to the gradient of p p . 
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FIG. 3: (Color online) The energy per flux quantum E n /n, in units of Eb os (see Eq. (19)), as a function of the number n of units 
of flux in the flux tube. Top left, simple proton superconductor with neutrons completely decoupled (/3 = a = 0); top right, 
density coupling between condensates (/3 = .5, a = 0); bottom left, gradient coupling between condensates (/? = 0, a = .5); 
bottom right, both couplings (/? = a — .5). 



B. Energetic stability of flux tubes 



In Fig. 3, the energy per flux unit {E n /n) is plotted against n for various values of the Ginzburg-Landau parameters, 
namely re, the density coupling f3, and the gradient coupling a. We fixed (</>p) 2 /(</>n) 2 = 0.05 (Sec. IIIA). 



1. No coupling to neutrons 

The upper left plot of Fig. 3 shows En C ^ (re) / n, the energy per flux quantum when there are no interactions between 
the neutron and proton pairs. We see that the only possible phases are the standard type I and type II, with a 
transition at the Bogomolnyi point, re = 1/v2j where the favored value of n jumps from 1 to infinity.. The lower line 
(re just below l/y/2) corresponds to type-I, where the lowest energy/flux is at n = oo, so flux tubes attract. The upper 
line (re just above l/y/2) corresponds to type-II, where the lowest energy/flux is at n = 1, so flux tubes always repel 
each other. The middle line corresponds to the transition point (re = l/y/2), where there is no interaction between 
flux tubes [25]. Our numerical results are consistent with the expected form (19): when <5re > the asymptotic value 
of E n /n is increased, and E n /n rises monotonically towards that asymptotic value, and conversely when Sn < the 
asymptotic value of E n /n is decreased, and E n /n falls monotonically towards that asymptotic value. It is clear that 
ci in (19) must be positive to obtain this behavior at large n. From fits to our numerical calculations we find that c\ 
is always positive, so it "fights against" the leading ci/y/n term, but for all n }z 1 it is overwhelmed. In fact, we find 
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that (19) gives an excellent fit to our results down to n — 1, without any higher order terms. 

In the remaining panels of Fig. 3, we explore the effect of density and gradient couplings between the proton 
superconductor and the neutron superfluid. 



2. Density coupling to neutrons 



The upper right panel of Fig. 3 shows the effect of a density coupling between the condensates. From (19) and (20) 
we expect 

E n /n = E Boe + (MSk - Mp) + 2 2 - + + • • • (22) 

The first point to notice is that the density coupling shifts the critical k to a larger value. The transition between type-I 
and type-II occurs when the asymptotic behavior at large n changes from rising to falling, i.e. when the coefficient of 
the 1 / y/n term changes sign. This occurs for some positive value of 5 k 

M g bi 

(fccritGS) = "Tp 1 cc P 2 (23) 
Mci 

2 

which rises as (3 2 because M, Mp, bi, and ci arc all positive, and Mp oc /3 2 when a = (Sec. IV A 4). Thus in the 
upper right panel of Fig. 3 we had to increase k from around 0.707 to around 0.818 in order to find the transition. 

The other important point is the presence of a minimum in E n /n when k is just above the new type-I/type-II 
boundary, indicating that the favored value of n may be neither 1 (standard type-II) nor infinity (type-I) but some 
intermediate value. This is consistent with (22), as long as we assume that the coefficient b\ from (20) is either positive, 
or negative and of sufficiently small magnitude, so that the 1/n term in (22) has a positive coefficient (recall that 
Mp, M, and c\ are all positive, and 5k is also positive in this region). The minimum will then arise from competition 
between the positive 1/n term, which dominates at smaller n, giving a negative slope, and the 1/y/n term which has a 
negative coefficient (because 5k is just above the new critical value) and dominates at larger n giving a positive slope. 
However, as 5k is reduced the negative coefficient of 1/y/n becomes smaller and smaller, and the minimum moves out 
to arbitrarily large n, so the energetically favored value of n does not jump suddenly from 1 to oo as in the standard 
case, but increases in steps from 1 to infinity as we lower k through a range of values down to the new critical value. 
This creates an infinite number of "type-II(n)" phases, each with a different flux in the favored flux tube, and when 
that flux becomes infinite the superconductor becomes typc-I. This behavior is seen in our numerical results (Fig. 4). 



3. Gradient coupling to neutrons 

The lower left panel of Fig. 3 shows the effect of a gradient interaction with the superfluid. From (19) and (20) we 
expect 

—M„Sl — SkMci —M Si 4- 5kMc-\ 

E n /n = E Bos + M5K+ 1-+ Mg«i + o*Mc 1 + _ (24) 

y/n n 

Here we see that the gradient coupling shifts the critical k to a smaller value. The transition between type-I and 
type-II occurs when the coefficient of the 1/y/n term changes sign, which in this case happens for small negative 5k, 

M a si 

^Kcrit(o') = — — —2- oc -a (25) 
Mci 

2 

which is proportional to —a 2 because M, M a , si, and ci are all positive, and M a cx a 2 when j3 = (Sec. IVA4). 

The other important feature of this plot is the presence of a maximum in E n /n when k is close to the type-I/type-II 
boundary. This is consistent with (24), as long as we assume that the coefficient Si from (21) is either positive, or 
negative and of sufficiently small magnitude, so that the 1/n term in (24) has a negative coefficient. The maximum 
will then arise from competition between the negative 1/n term, which dominates at smaller n, giving a positive slope, 
and the 1/y/n term, which dominates at larger n giving a negative slope. 

The presence of this maximum allows for the possibility of metastable flux configurations. If we scan down in k, we 
start in a type-II region where E n /n has its minimum at n — 1 and rises monotonically with n. But at some point 
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a metastable minimum at n = oo appears, which drops to become degenerate with the minimum at n = 1. At this 
point there is a first-order transition: at the critical field, n = 1 flux tubes would co-exist with macroscopic normal 
regions (i.e. flux tubes with n = oo) but not with flux tubes of intermediate size. Reducing k further, the n = 1 flux 
tube becomes energetically metastable, and finally unstable. 

4- Density and gradient coupling to neutrons 

The lower right panel of Fig. 3 shows the effect of a combination of gradient and density interactions. As k is 
decreased, a metastable energy minimum emerges at finite n; it drops and becomes a new global minimum at n = n* , 
yielding a sharp transition from n = 1 typc-II to n — n* type-II. As k is reduced further the favored number of flux 
quanta in a flux tube rises in integer steps from n* to infinity, at which point the superconductor becomes type-I. 



Density coupling to neutrons 



Density coupling to neutrons 
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FIG. 4: (Color online) Effect on the superconductor of density coupling j3 to a superfluid, displayed as a phase diagram in the 
K-/3 plane, with no gradient coupling {a = 0) and (4>p) 2 / {<pn} 2 = 0.05. The left panel shows how non-zero (3 causes an increase 
in /tcritical- In the right panel we magnify the transition region near j3 = 0.5, illustrating that on the type-II side there is a 
sequence of "type-II(n)" bands in which the number of flux quanta in the favored flux tube rises, reaching infinity when the 
superconductor becomes type I. 



Gradient coupling to neutrons 



Type-II (meta-I) 




FIG. 5: (Color online) Effect on the superconductor of gradient coupling a to a superfluid, displayed as a phase diagram in the 
K-a plane, with no density coupling (j3 — 0) and (</ > p) 2 /(</'n) 2 = 0.05. The gradient coupling causes a decrease in « C ritical, and 
creates metastable states on either side of the transition, with spinodal lines as shown. 
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FIG. 6: (Color online) Phase diagram for combined density and gradient interactions: the K-0 plane for a — 0.5 and 
(4> p } 2 1 '((j> n ) 2 = 0.05. The type-I/type-II boundary is no longer symmetric under f3 — ► — (3. In the right panel we magnify 
the transition region near (3 = 0.5, illustrating that on the type-II side as k, decreases the number of flux quanta in the favored 
flux tube jumps from 1 to a finite value (in this case n = 5) and then there is a sequence of bands in which n rises to infinity, 
at which point the superconductor becomes type I. 




FIG. 7: (Color online) Phase diagrams in the k vs. (<f> n ) 2 /(4> p ) 2 plane. Vertical dashed lines show {4>n) 2 / {<j> v ) 2 = 20, the value 
used for other figures in this paper. The left panel is for density coupling (3 = —0.1, but no gradient coupling (a = 0). The 
right panel is for gradient coupling a — 0.1, but no density coupling (/3 = 0). In both cases, we see that the type-I/type-II 
transition converges to ft = 1/V% as the neutron condensate disappears. For the case of a density coupling, as the neutron 
condensate decreases, the type-I/type-II boundary changes at (4>n) 2 / (4>p) 2 ~ 10 from a narrow region of type-II(n) phase bands 
(thick line) to wider metastable regions. 



C. Phase diagrams 



Figures 4-7 illustrate the additional structure in the phase diagram of the superconductor induced by the couplings 
to a superfluid. Each diagram is a two-dimensional slice through the parameter space. 

Figure 4 shows the consequences of a density coupling (3 between the superfluid and superconductor. We see that 
the density coupling, irrespective of its sign, favors typc-I superconductivity, pushing the the critical n for the type- 
I/type-II transition up to higher values, forming a parabolic phase boundary in the (3-k plane, as expected from (23). 
This can be thought of as arising from the fact that nonzero (3 lowers the energy per flux of the core of large flux 
tubes (see (22)), which favors type-I superconductivity. 

In the right panel we zoom in on the transition line near (3 = 0.5 to show the substructure in the phase transition 
region that is invisibly small in the left panel. As one would expect from our discussion of Figure 3 (upper right 
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panel) , on the type-II side of the transition there is a series of bands distinguished by the number of flux quanta n in 
the energetically favored flux tube. "Type-II (n = 1)" is the standard type-II superconductor. With decreasing re wc 
find transitions to Type-II (n = 2), Type-II (n = 3), and on up to n = oo which is a typc-I superconductor. 

In Figure 5 we show the consequences of a gradient coupling a between the supcrfluid and superconductor. We 
see that the gradient coupling, irrespective of its sign, favors type-II superconductivity, pushing the critical re for the 
type-I/type-II transition down to lower values, forming an inverted parabolic phase boundary in the er-re plane, as 
expected from (25). It also makes the phase transition first order, with spinodal lines where the unfavored phase 
becomes metastable. Both these effects arise from the lowering of the energy of the wall of the vortex, as explained 
in Sec. IV B 3. 

In Figure 6 we show phase diagrams for the combination of both density and gradient couplings, fixing a = 0.5 and 
varying [3. As discussed in Sec. IV A 4, we expect that when a ^ the [3 — > —(3 symmetry is now broken. In the right 
panel we magnify the transition region near /3 = 0.5, illustrating that on the type-II side as re decreases the number 
of flux quanta in the favored flux tube jumps from 1 to a finite value n = 5, and then there is a sequence of bands in 
which n rises, reaching infinity when the superconductor becomes type I. This is the expected behavior, based on our 
discussion in Sec. IV B 4. 

Finally, in Figure 7, we anticipate one direction in which this work could be extended, by exploring the consequences 
of varying the ratio of the superfluid density to the superconductor density, which up to now was fixed to (4> n ) 2 / (4>p) 2 — 
20, an appropriate value for neutral beta-equilibrated nuclear matter, of the type we expect to find inside neutron 
stars. Figure 7 shows phase diagrams in the plane of re and (<fi n } 2 / {<j) p } 2 for a system with a density coupling (left 
panel) and with a gradient coupling (right panel). 

For the case of a density coupling we use a negative value of the coupling, because this corresponds to an attractive 
interaction, which gives smooth behavior in the limit where the neutron condensate disappears, (4>n} 2 / (4>p) 2 ~ * 0. 
As is clear from the plot, the type-I/type-II transition then converges to the standard value for a single-component 
superconductor, re = \/\/2. For a repulsive interaction, the ((f> n ) 2 / (4>p) 2 ~~* limit is singular: we discuss this in more 
detail below. It is interesting to note that the effects of the density coupling change dramatically with the relative 
densities of the neutrons and protons. At {<f> n ) 2 / ' (4> p ) 2 > 10 the density coupling produces a thin region of multi-flux- 
quantum "type-II(n)" phases, as was illustrated in Fig. 4. But for lower values, it has a similar effect to a gradient 
coupling, inducing metastable regions on either side of the type-I/type-II boundary. This should be understandable 
in terms of the dependence of the coefficients bi and b\ (Eqn. (20)) on (<fr n ) 2 / (<fr p ) 2 . In Sec. IVB2 we argued that 
if bi is large enough then the E n /n curve has a minimum at finite n, yielding a typc-II(n) phase. We conjecture 
that as (4> n ) 2 / (4>p) 2 gets smaller, b\ becomes sufficiently negative that this is no longer the case, and instead there is 
a maximum, leading to mctastability of the n = and n — oo states in spinodal regions around the type-I/typc-II 
boundary. This is a topic for future investigation. 

For the case of a gradient coupling (right panel of Fig. 7) the effects of varying ((f> n ) 2 / (4> P ) 2 are less dramatic. It 
is interesting that, as for a density coupling, the variation is non-monotonic. Again, we conjecture that this could be 
understood in terms of variation of the coefficients si and si (Eqn. (21)) with (<fi n ) 2 / (4> p ) 2 . As the superfluid density 

drops to zero, its effects become negligible, and the critical value of re converges towards l/y/2 as one would expect. 

Finally, we discuss the singularity of the (4>n} 2 / (4>p) 2 ~ * limit for a positive (3, i.e. a repulsive density coupling 
between the neutron and proton condensates. From (6) we see that the expectation value of the neutron condensate 
is (</> n ) + \f3({4> p ) — 4> p ), so far from the flux tube, where <p p is (4> p ), it is (<fin)- But in the core of the condensate it is 
larger (there is less proton condensate to repel it). In fact, even if the parameter (4> n ) 2 were zero or slightly negative, 
there would be a positive neutron condensate in the core of the flux tube. This shows that for positive j3 the neutrons 
do not decouple and become irrelevant in the limit (</>„) — » 0. We note two consequences of this. Firstly, for small 
(4> n ) the (3 — > — (3 symmetry discussed in Sec. IV A 4 is no longer present, because the effect of the flux tube on the 
neutron condensate is no longer a small perturbation. Secondly, in a system where {4> n ) 2 is small and negative (i.e. the 
neutrons just barely fail to condense in the presence of the proton condensate) flux tubes could have supcrfluid cores, 
which is another topic that we leave for future investigation. 

V. CONCLUSION 

We conclude that coupling a superconductor to a co-existing supcrfluid causes significant modification of the ener- 
getics of the flux tubes. On the basis of calculations restricted to the cylindrical geometry of n-quantum flux tubes, 
we conclude that a coupling between the densities of the condensates shifts the type-I/type-II boundary to larger re, 
and, if the superfluid density is high enough, appears to create an infinite number of new "type-II(n)" phases whose 
most stable flux tubes contain multiples of the basic flux quantum. A gradient coupling between the condensates 
leads to metastable regions surrounding the transition between type-I and type-II superconductivity. 
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As discussed in Section II, our calculation corresponds to comparing the energy at zero and infinite separation 
of flux tubes with varying numbers of flux quanta. This leaves open the possibility that there might be additional 
minima at finite separation. It is therefore possible that in parts of the phase diagram there might be a different 
phase from the ones we identify, namely an alternative type of type-II superconductor in which the spacing between 
flux tubes is fixed by the microscopic physics rather than by the strength of the applied field. To resolve this question 
will require calculation of the free energy of a pair of flux tubes at arbitrary separation. Such calculations have 
been performed for large separation [28-30], and by perturbing about the Bogomolnyi point [17] and by numerical 
computation [31]. In particular, the numerical methods that have been used recently to follow the interaction and 
annihilation or vortcx-antivortex pairs [32] would be readily applicable to the simpler time-independent calculation 
of the interaction potential of flux tubes. Another natural generalization of our calculation would be to allow for 
non-s-wave pairing, such as the 3 P2 pairing that is believed to occur in the neutron superfluid in the core of a neutron 
star. 

Our results add another example to the class of two-component Ginzburg-Landau models with non-standard su- 
perconducting behavior. Previous work in this area includes the SO(5) model of high-temperature superconductivity, 
which has flux tubes described by a two-component GL model, where each component carries a different U(l) charge, 
and only one of them condenses in the vacuum [13]. Another example is the case of a two-component GL model where 
both components have electric charge, very different mass, and nearly the same Fermi energy. This system was found 
to have non-monotonic E(n)/n and intermediate minima in the interaction potential [7]. 

The exotic phenomena that we predict are localized to the region around the type-I/type-II transition, so they may 
not turn out to be relevant for the inner core of a neutron star, which is believed to be well inside the type-II regime 
[1]. However, given the extremely impressive recent progress in creating exotic systems such as multi-component 
superfluids of trapped cold atoms, it seems quite conceivable that a material that is both a superconductor and a 
superfluid might be created in the laboratory, and could be studied under controlled conditions. Our results would 
be directly relevant to such a material. 
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